- Descriptive point pattern analysis
- Kernel density maps
- Polygon extraction
We are going to use the tree species dataset from two weeks ago, more specifically the occurences of abies alba (silver fir):
library(sp)
library(raster)
abies <- read.csv("data/abiesalba.csv")
head(abies)
## long lat ## 1 -0.01377021 48.78623 ## 2 -0.01423346 42.83379 ## 3 -0.02027734 42.79684 ## 4 -0.02039473 43.01494 ## 5 -0.02171301 48.48641 ## 6 -0.02336757 42.92377
abies_sp <- SpatialPoints(abies, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
abies_sp <- spTransform(abies_sp, CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
The centroid can easily be calculated by averaging the coordinates:
coords <- abies_sp@coords coords <- as.data.frame(coords) xmean <- mean(coords[, 1]) ymean <- mean(coords[, 2]) p <- ggplot(coords, aes(x = long, y = lat)) + geom_point(alpha = 0.3) + geom_point(aes(x = xmean, y = ymean), col = "red", shape = 2)
We can also calculate the median of the coordinates:
xmedian <- median(coords[, 1]) ymedian <- median(coords[, 2]) p <- ggplot(coords, aes(x = long, y = lat)) + geom_point(alpha = 0.3) + geom_point(aes(x = xmean, y = ymean), col = "red", shape = 2) + geom_point(aes(x = xmedian, y = ymedian), col = "blue", shape = 2)
We can calculate the average/minimum/maximum distance between points:
dist <- raster::pointDistance(abies_sp, longlat = FALSE) dist[1:3, 1:3]
## [,1] [,2] [,3] ## [1,] 0.0 661483.862 665579.331 ## [2,] 661483.9 0.000 4123.066 ## [3,] 665579.3 4123.066 0.000
dist <- dist[lower.tri(dist)] c(min(dist), mean(dist), max(dist)) / 10000
## [1] 0.04996316 55.87651901 268.64010711
We can generate a grid from the points extent::
ext <- extent(abies_sp) ext
## class : Extent ## xmin : 3185500 ## xmax : 5621500 ## ymin : 1690500 ## ymax : 3669500
ras <- raster(ext) ras
## class : RasterLayer ## dimensions : 10, 10, 100 (nrow, ncol, ncell) ## resolution : 243600, 197900 (x, y) ## extent : 3185500, 5621500, 1690500, 3669500 (xmin, xmax, ymin, ymax) ## coord. ref. : NA
And set the resolution to 1000*1000 meters while keeping the projection:
dim(ras) <- c(1000, 1000, 1) ras
## class : RasterLayer ## dimensions : 1000, 1000, 1e+06 (nrow, ncol, ncell) ## resolution : 2436, 1979 (x, y) ## extent : 3185500, 5621500, 1690500, 3669500 (xmin, xmax, ymin, ymax) ## coord. ref. : NA
projection(ras) <- projection(abies_sp)
Following we can use this raster to rasterize the points and count the number of points in each cell:
countgrid <- rasterize(abies_sp, ras, fun = "count", background = 0) mapview(countgrid)
We can easily calculate the distances to cells with Abies. However, I first aggregate the raster to a loer resolution (otherwise it will take too long):
abies_ras <- aggregate(countgrid, fact = 10, fun = sum) mapview(abies_ras)
Following I set all cells containing no Abies to NA and calculate the distance:
abies_ras[abies_ras == 0] <- NA distance_abies <- distance(abies_ras) mapview(distance_abies)
Kernel density estimates do not simply count the occurences, but 'smooth' by fitting a function (e.g., multivariate-normal) to the data. We can calculate Kernel densities as raster:
library(spatialEco) kernel_ras <- sp.kde(abies_sp, bw = 10000, newdata = ras) mapview(kernel_ras)
ggplot2 can also estimate kernels and nicely plot them atop of Google Maps!
library(ggmap) map <- get_map(location = c(lon = mean(abies$long), lat = mean(abies$lat)), zoom = 5)
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=47.077251,9.542565&zoom=5&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
p <- ggmap(map) +
geom_density2d(data = abies, aes(x = long, y = lat), size = 0.3) +
stat_density2d(data = abies,
aes(x = long, y = lat, fill = ..level.., alpha = ..level..), size = 0.01,
bins = 50, geom = "polygon") +
scale_fill_gradient(low = "green", high = "red") +
scale_alpha(range = c(0.1, 0.5), guide = FALSE)